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; Abstract 
5— i ■ 

Recently, a new approach for the stabilization of the incompressible Navier-Stokes equations for higher 
Reynolds numbers was introduced based on the nonlinear differential filtering of solutions on every time step 
of a discrete scheme. In this paper, the stabilization is shown to be equivalent to a certain eddy-viscosity 
C^) . model in LES. This allows a refined analysis and further understanding of desired filter properties. We also 

consider the application of the filtering in a projection (pressure correction) method, the standard splitting 
algorithm for time integration of the incompressible fluid equations. The paper proves an estimate on the 
convergence of the filtered numerical solution to the corresponding DNS solution. 

. 1 Introduction 

A stabilization of a numerical time-integration algorithm for the incompressible Navier-Stokes equations 

u t + (u- V)u - vAu + Vp = / in „ x ( T] (1) 

_ , div u = 

CN 

\ for large Reynolds numbers with the help of an additional filtering step was recently introduced in pQ . Denote by 
00 ■ w n or u n approximations to the Navier-Stokes system velocity solution at time t n , and similarly p n approximates 
pressure p(t n ). Let At = t n +i — t n . The algorithm, referred to further as (Al), reads: For n = 0, 1, . . . and 

■ u° = u(t°) 
(N . 

, 1. compute intermediate velocity w n+ from 

m 

-(w n+1 - u n ) + (w n+1 ■ V)w n+1 + Wp n+1 - vAw n+1 = f n+ \ 



At 

subject to appropriate boundary conditions; 



divw n+L = 0, 



2. filter the intermediate velocity, w n+1 := Fw n+1 ; 



3. relax u n+1 := (1 — x) wn+1 + X w7l+1 i with a relaxation parameter \ G [0, 1]. 

Here F is a generic nonlinear filter acting from L 2 (£l) 3 to i7 1 (fi) 3 . We shall consider further in the paper 
several examples of differential filters. The convergence of the finite element solutions of (Al) to the smooth 
Navier-Stokes solution has been analyzed in pQ. One advantage of the approach is the convenience of imple- 
mentation within an existing CFD code for laminar flows and flexibility in the choice of a filter. Numerical 
results from [21 [SI HI SI IS] with composite nonlinear differential filters, as defined in Section [31 consistently show 
more precise localization of model viscosity and its more precise correlation with the action of nonlinearity on 
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the smallest resolved scales than plain Smagorinsky type LES or VMS methods. Thus we deem the approach 
deserves further study, should be put into perspective and related to developing LES models. 

In this paper, we show that introducing the filter stabilization is closely related (and even equivalent in 
a sense which is made precise further in the paper) to adapting a certain eddy-viscosity model for LES. The 
connection to a LES model helps us to quantify the model dissipation introduced by the filter stabilization 
(Theorem [1]), formulate stability criteria (see ([6]) and (jH])), and gives insight into the choice of the filter and 
the relaxation parameter. In particular, it provides an explanation why the stabilization by the filtering avoids 
adding excessive model viscosity in regions of larger velocity gradients, unlike most other eddy viscosity models. 

The entire approach is specifically designed for treating higher Reynolds number flows. Therefore, it is 
natural to extend it to the Chorin-Temam-Yancnko type splitting algorithms, which are the prevailing method 
for the time-integration of the incompressible Navicr-Stokcs equations for fast unsteady flows. Such (rather 
natural) extension is presented in the paper together with the relevant error analysis. We note right away that 
the analysis demonstrates the convergence of numerical solutions to the Navier-Stokes smooth solution, while 
it would be also interesting to analyze the error of the numerical solutions to a (presumably smoother) solution 
of the corresponding LES model. However, the specific difficulty we faced in the latter case is the lacking of the 
monotone property by most of eddy viscosity indicator functionals, which were numerically proved to be useful 
in defining the filter F, see Section [3] Though practically attractive, introducing such functionals makes the 
mathematical well-posedness of the LES model and accordingly the error analysis hard to accomplish and we 
are unaware of relevant results in this direction. 



2 Filter stabilization and LES model 

It is well known, see, e.g., [6] or [7], that explicit filtering is related to adding eddy or artificial viscosity. The 
connection of the filter stabilization as defined above to LES modeling is easily recovered by noting that shifting 
the index n + 1 — > n on steps 2 and 3 and using step 1 gives the implicit discretization of the Navier-Stokes 
equations, with explicitly treated nonlinear dissipation term: 

-^-(w" +1 - w n ) + {w n+1 ■ V)w n+1 + Vp n+1 - vAw n+1 + -^Gw n = f n+1 , 

divw n+1 = 0, 
with 

G := I — F, I is the identity operator. 

Assume x — XoAi, where xo is a time- and mesh-independent constant, then ((2|) can be treated as the time- 
stepping scheme for 

w t + (w ■ \7)w + Vp- uAw + xoGw = /, 

div w = 0. 

These arguments show that the numerical integrator (Al) with filter stabilization is the splitting scheme for 
solving ((3|). Furthermore, ((3]) can be observed as a LES model, with xoGw corresponding to the Reynolds stress 
tensor closure: 

V ■ (w g>w — w (g)w) ~ xoG w. 

This simple observation leads to a refined analysis and better interpretation of the numerical results and the 
method properties. 

We note that x = O(At) is exactly the scaling of relaxation parameter which allows us to prove optimal 
convergence result for a time-stepping splitting method (Theorem [3|). Furthermore, numerical experiments in 
3, 8 suggested that x = O(At) is indeed the right scaling of the relaxation parameter with respect to numerical 
solution accuracy. 

We start by showing several numerical properties of the approach. Throughout the paper we use (•, •) and j| • || 
to denote the L 2 scalar product and the norm, respectively. For the sake of analysis, assume the homogeneous 
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Dirichlct boundary conditions for velocity. Taking the L 2 scalar product of (J2]) with 2Atw n+1 and integrating 
by parts gives 

|| w «+l||2 _ || w n||2 + l\\ w n+l _ w n\\2 + l/ At\\Vw n+1 \\ 2 + X (Gw 7 \ W n+1 ) = At(f 1+1 , W n+1 ) . (4) 

For a self-adjoint filtering operator, i.e. (Gu,v) = (Gv,u) for any u,v £ Hq(£1) 3 , the equality §4§ can be 
alternatively written as 

|K +1 || 2 - |K|| 2 + j/Ai||Vu;' l+1 || 2 + | ((Gw n+ \w n+1 ) + (Gw 7 \w n )) 

= At(f, w n+1 ) + i ( x (G(w n+1 - w n ), w n+1 - w n ) - \\w n+1 - w n f) . (5) 

Considering the last two terms on the right-hand side, we immediately get the sufficient condition of the energy 
stability of © for the case of self-adjoint filters: 

X (Gu,u) < |M| 2 VueffoW- (6) 
If G is not necessarily self-adjoint, one may rewrite (|4|) as 

||w" +1 || 2 - ||w"|| 2 + |ll w ' rl+1 -t^ ?l H 2 + ^^*l|Vw; n+1 || 2 + x(G'w rl ,w") = At(/, + x(Gw",-u; Tl - 

Thanks to the Cauchy inequality one gets for any 6 £ M: 

|K +1 || 2 - ||w"|| 2 + j/At||Vu- n+1 || 2 + (1 - e) X (Gw n ,w n ) 

< At(f, w n+1 ) - X (e(Gw n ,w n ) - ^(Gw n ,Gw n )) . (7) 

In this more general case, one may consider the following sufficient condition for the energy stability. Fixing, 
for example, 9 = i, assures the sum of the last two terms in ([7]) is positive if 

X (Gu, Gu) < (Gu, u) VueH^nf. (8) 

Assume G is self-adjoint and w n approximates a smooth in time Navier-Stokes solution, then (|5|) leads to 
the following energy balance relation of the numerical method: 

N N N 

\\w N \\ 2 + vY, AtHV^l 2 + X o£ At(Gw n ,w n ) = \\w°\\ 2 + £ At(f\ w n ) + O(At). 

n—l n—1 n—1 

In particular, we may conclude that the hltcr stabilization introduces the model dissipation of 

N 

Xo^A*(G W ", W "). (9) 

n=l 

Finally, we notice that the filtering and relaxation steps in (Al) can be rearranged as 



-XoGw 



n+l 



At 

which is the explicit Euler method for integrating 

u t = -XoGu on [t n ,t n+1 ], with u(t n ) = w(t n+1 ). (10) 



The coupling of a DNS method with the evolution equation (|10|) is known as another way of introducing explicit 
filtering in modelling of dynamical systems, e.g. [6]. This suggests that an improvement leading to higher order 
methods for integrating (|10[) might be possible. 

In the next section, we shall study properties of the operator G for a class of nonlinear differential filters. 
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3 Nonlinear Differential Filters 



Linear differential filters have a long history in LES, see [9]. We also point to [TO] and references therein for 
applications of linear differential filters in the Lagrange-averaging turbulence models. In this section, we consider 
a family of nonlinear differential filters for the filtering procedure. Some conclusions will be drawn concerning 
the stability conditions ((6]), ((8j) and equivalence to other approaches in the LES modelling. We use the following 
notation: 

V := {v e H^{n) 3 : divv = 0} , H = {v e L 2 (fL) 3 : divw = 0,v ■ n\ dn = 0} . 

By P wc denote the L 2 orthogonal projector from L 2 (il) 3 onto H. 

For a given sufficiently smooth vector function u and w € L 2 (f2) 3 we define Fib as the solution to 

{S 2 a{u) V(F w), Vv) + (F w, v) = (w, v) \fv G X, (II) 

with an indicator functional < a(u) < 1 and filtering radius S 2 , which generally may depend on x and t, 
<5max = maxj^ \S\. Here X = H^ft) 3 or X = V, if the filter is div-free preserving. We note that it is not 
immediately clear if the problem (|II[) is well-posed. In practice, this is not an issue, since in a finite dimension 
setting, e.g. for a finite element method, the bilinear form from the left-hand side of (fTTj) is elliptic and thus 
(fTTj) is well-posed. Otherwise, we may assume < e < a(u) < 1 for some sufficiently small positive e. If we 
assume this, none of our results further in the paper depend on the parameter e. It is standard to base the 
indicator functional on the input function w itself, that is u = w and we will denote W := F w in this case. 
However, in the course of analysis we need to consider (auxiliary) filtering with u ^ w. If we need to show 
explicitly the function used for the indicator, we shall write F(u)w instead of Fw or F(w)w instead of w. 
The action of G = I — F, w g := G w, is defined formally as the solution to 

(5 2 a(u)\7w g , Vv) + {w g ,v) = (5 2 a{u)Vw, Vv) W £ X. (12) 

The operator G is self-adjoint on X and in the operator notation it can be written as 

G = -[I- AJ-'A,, (13) 

with 




div(<5 2 a(u)V) if X = H^{n) 3 , 
Pdiv((5 2 a(w)V) if X = V. 



Since operator A Q is self-adjoint and positive definite, one see from ([T31 that G < I and thus the sufficient 
stability condition (|6|) holds for any x G [0,1]- This can be easily verified in a formal way by substituting v — Fw 
in dTTJ) to get (w,Fw) > and thus (w,Gw) = (w,w — F w) < \\w\\ 2 for any w £ Hq(H) 3 . Moreover, varying 9 
in and using ((HJ, one shows the energy stability estimate for any \ £ [0, 2]. However, such refinement is not 
important for our further analysis. 

With the help of © and (|13p . we now quantify the model dissipation introduced by the differential filters. 
To make notation shorter and without loss of generality, let \ = Xo^^- 

First, representation (|13[) immediately implies G < — A a . Thus the additional dissipation introduced by the 
differential filtering does not exceed those introduced by the LES closure model: 

div(w Cg> w — w (g> w) w — xo A a u>. (14) 

It is easy to show that for a discrete case and if the condition 

5 < spatial mesh width 

holds and < a(u) < 1, then the dissipation introduced by the differential filtering (|lll) is equivalent to the 
dissipation of the closure model (JUJ). 

We make the above statement more precise for a finite element discretization. To this end, assume a 
consistent triangulation T of CI, satisfying the minimal angle condition 

inf p{K)/r{K) =: a > 
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where p(K) and r(K) are the diameters of inscribed and superscribed circles (spheres in 3D) for a triangle 
(tetrahedron) K. We have the following result. 

Theorem 1 Assume X is the finite element space of continuous functions which are polynomials of degree 
p > 1 on every element K and ni&x X £K \S(x)\ < C$ r(K) for any K El~, with a constant Cs independent of K . 
Then for any w G X the equivalence 

c(S 2 a(u)Vw,\7w) < (Gw,w) < (5 2 a(u)Vw,\7w) (15) 

holds with a constant c> independent ofw, the indicator a(-) , and the filtering radius 6. The constant c> 
may depend on p, Cs, and ao- 

Proof. Consider the finite element inverse inequality 

\\Vw\\l*(k) ^copWIHU'Ck) VweX, (16) 

where the constant Co depends only on the polynomial degree p and ao- The inequality (fT6"|) . the assumption 
on <5 and the minimal angle condition imply 

\\^w\\ L 2 (K) < C\\w\\ L 2 {K) , (17) 

where the constant C depends only on p, Cs, and ao- Squaring ([T7|). summing over all K G T , and recalling 
that a(-) < 1, implies 

(5 2 a(u)Vw,Vw) < C 2 \\w\\ 2 . (18) 
Denote w g = Gw for some w G X. We set v = w g and v = —w in f|12[) and sum up the equalities to get 

= (S 2 a(u)Vw g , Vw g ) + (w g ,w g ) - 2 ((J 2 a (ti) Vik, Vw 9 ) - (w g , w) + (S 2 a(u)Vw, Vw) 
= \\w g \\ 2 — {w g ,w) + (5 2 a(u)S7(w - w g ),V(w - w g )). 

Thus, it holds ||w ff || 2 < (w g ,w), i.e. the condition ((HI). Now we set v = w in (fT2|) and use ([8]) and (fT8j) to 

estimate 

(S 2 a(u)Vw, Vio) = (5 2 a(u)Vw g , Vw) + (w g , w) 

< ^(S 2 a(u)\7w g , Vw g ) + ^(S 2 a(u)\7w, Vw) + {w g ,w) 

< lc 2 \\w g \\ 2 + ±(5 2 a(u)\7w,Vw) + (w g ,w) 

< i\c 2 + l){w g ,w) + ^{5 2 a(u)\7w,Vw). 

We proved the lower bound in ([TSj) . 

To show the upper bound we set v = w g and v = w in (|12p and sum up the equalities to get 

= (S 2 a(u)\7wg, Vw g ) + (w g ,w g ) + (w g , w) — (^a^Vw, Vw). 

This yields the upper bound in (fT5|) : (w g ,w) < (<J 2 a(«)Vw, Vi«). 

Few conclusions can be drawn from the equivalence result (|15[) concerning the relation of the filter stabiliza- 
tion to some other eddy-viscosity models. 

The use of the linear differential filter (a = 1), as considered in [3], is equivalent to the method of artificial 
viscosity. This means that the model dissipation is equivalent to the isotropic diffusion scaled with xo$ 2 - Given 
what is known about the method of artificial viscosity, it is not surprising that the method is not very accurate 
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in this case. Thus, more elaborated indicator functionals should be used. Generally, we may think of a(u) as a 
real valued functional, depending on u, Vu, and selected with the intent that 



a(u(x)) ps for laminar regions or persistent flow structures, 
a(u(x)) w 1 for flow structures which decay rapidly. 



The choice of the Smagorinsky type indicator function, a(u) = |Vw|, does not necessarily satisfy the condition 
a(u) < 1. In this case, we do not have the equivalence result of the filter stabilization to the Smagorinsky LES 
model. Only the upper bound in (|15[) is guaranteed to hold. Thus the dissipation introduced by the filtering 
with a(u) = |Vu| is likely less than that of the Smagorinsky model. This can be a desirable property, since the 
Smagorinsky LES model is known to be severely over-diffusive for certain flows, e.g. [11] . and several ad hoc 
corrections were introduced such as van Driest damping, dynamic models, and others, see [12 ] 113 ] [14] . 

Several reasonable indicator functions a(u) are known to satisfy the boundedness condition: < a(u) < 1. 
These are the re- normalized Smagorinsky type indicator [15] , the indicator based on the Q-criteria [16j and the 
Vreman indicators |17j ; also an indicator based on the normalized helical density distribution was considered in 
[2]. Given several indicators cii(-), i = 1,. . . ,N, the combined indicator can be defined as the geometric mean: 



We remark, that the convergence results proved further in this paper do not rely on any smoothness prop- 
erties or particular form of a(-). 

The last remark in this section is that Theorem [T] does not give much insight if enforcing the divergence 
constraint in the filter is important or not. However, if we assume X — V in (|11[) . i.e., the filtered velocity 
satisfies the divergence free condition, then this slightly simplifies the error analysis in Section [5] 

4 Projection scheme with filter stabilization 

One idea behind introducing the filter stabilization or explicit filtering was to provide CFD software users 
and developers with a simple way to enhance existing codes for laminar incompressible flows to compute high 
Reynolds number flows. This goal is accomplished by making the filtering procedure algorithmically independent 
of a time integration method. Driven by this intention, we consider the Chorin [TS] splitting (projection) scheme 
with the additional separate filtering step. Projection methods are the common numerical approach to the 
incompressible Navier-Stokes equations and form a family of splitting algorithms, cf. [T9j[20j. We perform 
the numerical analysis for the simplest first order method given below. From the algorithmic standpoint, the 
generalization to higher order projection methods is straightforward, although analysis may become considerably 
more involved. 

Projection methods split the time evolution of the velocity vector field according to the momentum equation 
and the projection of the velocity to satisfy the divergence-free condition. The filtering step can be introduced 
before or after the projection step. In the former case, it is not necessary to augment the filter with the div-free 
constraint, since the projection step takes care of the keeping the approximates in the subspace of div-free 
functions. If the filter is div-free preserving, then it is reasonable to put it after the projection. In this paper 
we consider the constrained filter. We shall study the following algorithm: 

Step 1: Solve the convection-diffusion type problem: Given u n , w* , find w n+1 : 



The velocity w* is typically an interpolation from previous times, e.g. w* := w n or higher order interpo- 
lation. For the sake of analysis we consider w* = w n . 





■n+l 



(19) 



G 



Step 2: Project w n+1 on the div-free subspace: Find p n+1 and w n+1 solving the Neumann pressure Poisson 
problem: 

AV ' 

(20) 



+ Vp™ +1 


= 0, 


div w n+1 


= 0, 


w n+1 \ dn 


= 0. 



n+1 :=(l-x)w n+1 +Xw n+1 , (21) 



Step 3: Filter: w n+1 := Fw n+1 ; 
Step 4: Relax: 

with some x G [0) !]■ 

Similar to what was shown in section [21 shifting the index n + 1 — > n on steps 2-4 and substituting into (|19jl 
gives for \ = XoA< 

_ gS) + ( w * . v)w^+i + Vp n+1 - t/At^+i + xoG^ 1 - Atx GVp" +1 = f l+ \ 
Ai (22) 

diviyn+l - AtAp n+1 = 0. 

From (f2"2"j) we see that the splitting scheme (|19p - (|2"Tj) is formally the first order accurate time-discretization 
of the LES model ©. 

Further, we show that the splitting scheme (fT9j) - ((2Tj) is stable. There are two well-known approaches to 
accomplish the error analysis of projection methods. The one of Rannacher and Prohl [5D], [5T] uses the 
relation between projection and quasi-compressibility methods as it is seen from (|22[) . However, this analysis 
needs considerable effort to get extended to equations different from the plain Navier-Stokes equations. Another 
framework is mainly due to Shen (see [221 123] ). where convergence results were shown based on energy type 
estimates. In our error analysis we follow (to a certain extent) arguments from these two papers. 

5 Stability 

To show the stability of the splitting scheme, we need the following simple auxiliary result: 

Lemma 1 For w n+1 and u n+1 from the algorithm (fT!l]) - (|2ip and the filter F defined in (fTTj) . it holds 

IK +1 II > 

Proof. From the definition (fTTj) we obtain: 

1 



(d 2 a(iy n + 1 )Vw«+ 1 ,Vu; n + 1 ) + ||70 n + 1 || a = (w n+1 ,w" +1 ) = -(||u;" +1 || 2 + ||to n+1 || 2 - |K +1 - 1| 2 ). 
This yields 

|K +1 || 2 = 2(,5 2 a(u; ,l+1 )V^+ T ,V^ ;r + T ) + + Ww^ 1 - w n+1 f . 



Hence, ||-u/ i+1 || > From (J2TJ), we get 

\\u n+1 \\ < (i -x)IK +1 || + xlK+ T ll < IK +1 II for X £ [0,1]. 

Denote by || • the L 2 -dual norm for Hq(£1) 3 . Now we are ready to prove the following stability result. 
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Theorem 2 The algorithm (|19j) - (|2ip zs stable in the sense of the following a priori estimate: 
i-i i-i !-i i—i 

Ik'f + X! H w " +1 ~ ™" +1 l| 2 + XI ll^ 1 - u n \\ 2 + ^||Vu^+i|| 2 < ||w°|| 2 + ^ i/- 1 At||/(t„ +1 )|| 2 _ 1 (23) 

n— n— n— n— 

/or any i = 1, 2, . . . . 
Proof. 

Take the L 2 scalar product of (JTSJ with 2Atw™ +1 : 

2(u/»+ i - u",-^ 1 ) + 2j/At||V«^+' 1 || 2 = 2At(f n+1 ,^+ 1 ) < ^- 1 A<i|.f l+1 || 2 _ 1 + uAt\\ Vw^+i || 2 . 
Rewriting and simplifying this leads to: 

ll^ 1 !! 2 - IKH 2 + W^ 1 - u n f + j/At||V«^+i|| 2 < i/~ 1 Ai||/ n+1 ||i 1 . (24) 
The L 2 scalar of <[20j) with 2Atw n+1 and div w n+1 = gives 

2(w n+1 ~w^+\w n+1 ) = =► ||tu™ +1 || 2 -||i^|| 2 + |[iu" +1 -i^|| 2 = 0. 

Substituting ||u^+i|| 2 with ||w" +1 || 2 + - u^+i|| 2 in flM) yields 

|K +1 || 2 - ||u"|| 2 + -i^+i|| 2 + - ^f + i/AtHVto^H 2 < ^- 1 At||/' i+1 || 2 _ 1 . 

The application of Lemma [T] gives 

(|K +1 || 2 - ||u-"|| 2 ) + ||io"+ 1 - u^+i|| 2 + - u ll \\ 2 + vAtllVw^+if < v~ x A<||/ n+1 || 2 _ 1 . 

Summing up the inequality from n — 0, . . . , I — 1, we arrive at (|23l) . TJ 

6 Error Estimates 

We shall use (-, •) to denote the duality product between H~ s and Hq(Q) for all ,s > 0. In the following, we 
assume that the given data and solution to the equations ([TJ subject to the homogeneous Dirichlet velocity 
boundary conditions satisfy 

'u G {H 2 (9)) d C\V, 

f e L°°(0, T; (L 2 (n)) d ) n L 2 (0, T; (H\Q)) d ), 

f t eL 2 (0,T;H-% [ ' 

, su P*e[o,T] l|Vu(*)|| < C. 

We will use c and C as a generic positive constant which may depend on f2, v, T, constants from various Sobolev 
inequalities, uq, /, and the solution u through the constant C in (|25j) . 



Under the assumption (|25p one can prove the following inequalities, cf. 



sup {||u(t)|| 2 + ||«*(*)|| + ||Vp(t)||} < <7, (26) 

f£[0,T] 

i-T 

\\Vu t (t)\\ 2 +t\\u tt \\ 2 dt<C, (27) 

o 

which will be used in the sequel. Further we often use the following well-known [25] estimates for the bilinear 
form b(u, v, w) = J Q (u ■ V)v ■ w dx: 

c||V«||||Vi;||*||i;||ij||Viu||, 
b{u,v,w) < <( c||u|| 2 ||u||||Vu;]|, 
c||V«||||i;|| a ||«;||. 
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and b(u,v,w) = —b(u,w,v) for u G H. 

Define the Stokes operator Au = — FAu, Vu G D(A) = V PI H 2 (fl) 3 . We will use the following properties: 
A is an unbounded positive self- adjoint closed operator in H with domain D(A), and its inverse A -1 is compact 
in H and satisfies the following relations [22] [23] : 

f HA-^Ha < c||u|| and H^uH < c||«||v, 
3 c, C > 0, such that Vu G : < 

[ c\\u\\ 2 v ,<(A-\u)<C\\u\\ 2 vl . 

Before we proceed with the error analysis, we prove several auxiliary results given below in Lemma [2] The 
lemma gives estimates on the difference between a velocity w and the filtered velocity F(u)w. 

Lemma 2 Consider the differential filter F defined in (jlf I) with some sufficiently smooth vector function u. 
For any w G V and Fw G V it holds 

\\w-Fw\\ < WHViflH, (28) 

\\w-Fw\\y<6 2 max \\Vw\\. (29) 

Proof. Denote e = w — Fw. The equation gives 

{5 2 a(u) Ve, Vu) + (e, v) = (S 2 a(u)\7w, Vv) V v G V. 

Letting v = e yields 



||<5y / a(w)Ve|| 2 + ||e|| 2 = (5 2 a(u)Vw), Ve) < \\Sy/c^Vw\\\\5y/a(u)Ve\\ 

< l^v^HVell 2 + ±\\5^ajuj\7w\\ 2 < \\5^)Ve\\ 2 + \s 2 max \\Vw\\ 2 . 

This proves ((28)) . To show ([29]) , we note that setting v — F w — w in (|TT|) gives 

((J 2 o(u)VF)ii,V(fKi - w)) = - w\\ 2 < 0. 

Hence, we obtain: 

||(5V^t)VFu;|| 2 < \\Sy^ofu)^w\\ 2 . (30) 

Allowing v = A~ 1 {w — Fw) in (|TTj) leads to the following relations: 

||w-Fw||t>, = (w-Fw.^^w-Fw)) = {S 2 a{u)\7Fw,VA- 1 (w - Fw)) 

< ||(5 2 a(M)VFit)||||VA- 1 (u;-Fu;)|| < ^{\\5 2 a(u)WFw\\ 2 + \\w - Fw\\ 2 v ,) 

< ^Lxll^v/^yVF^H 2 + \\\w- Fw\\ 2 v ,. 

The last estimate and ([3T)|) implies ([2"9"]> . U 

Further in this section, we show that w n+1 , w n+1 and u n+1 are all strongly 0((Ai)s + <5) approximations 
to u(t n+ i) in L 2 (f2) 3 provided x = XoAi. Then we use this result to improve the error estimates to weakly 
0(At + S 2 ) approximations. This analysis largely follows the framework from [22] and [23] for the pure (non- 
filtered) Navier-Stokes equations, so we shall refer to these papers and [26] for some arguments which do not 
depend on the filtering procedure. 



Lemma 3 Let u be the solution to the Navier-Stokes system, satisfying (|25[) ■ Denote 

= u(t„ +1 ) - w^+\ e" +1 = u(i„ +1 ) - w n +\ and e n+1 = u(t n+1 ) - u n+1 . 
The following estimate holds 



l-i i-i 



WW 2 + E(H £ " +1 - e " +1 H 2 + - e "H 2 ) + E 2^At|| V6»+i|| 2 < C(At + C ax )- (31) 



n=0 



9 



Proof. Let R" denote the truncation error denned by 

^r(tt(t n +i) - «(*«)) - vAu{t n+1 ) + («(tn+i) ' V)«(t„+i) + Vjj(t„+i) = f n+1 + R n , (32) 
where R n is the integral residual of the Taylor series, i.e, 



1 



+1 



R n = ^- t J {t - t n )u tt {t)dt. 

By subtracting (|19p from (|32|) . we obtain 

i_ (e ^i _ e ») _ ^A^+i = (w n • V)^ 1 - («(t„ +1 ) ■ V)«(t„+i) - Vp(t„+i) + i?". (33) 
Taking the L 2 scalar product of ([33| with 2Ate n+1 , we get 

||c»+ 1 || 2 - ||e"|| 2 + He"^ 1 -e n || 2 +2i/At||Ve"+ 1 || 2 = 2A<(i? n , e^ 1 ) - 2Af(Vp(t„+i),e»+ 1 ) 

+ 2Atb*{w n ,^+ 1 ,e^ 1 ) -2Atb*{u{t n+1 ),u{t n+1 ),^). (34) 

The terms on the right-hand side are bounded exactly the same way as in p. 64 and [53] p. 512, leading to 
the estimates: 

At\b*(w n ,^,^)-b*(u(t n+1 ),u(t n+1 ),^)\ < ^\\V^\\ 2 +CAt\\e n f+C(At) 2 \\u t \\ 2 dt, (35) 

2At( J R",e™+ 1 ) < ^llV^f + C(At) 2 J* n+1 tlluttf^dt, (36) 

2At(Vp(Wi),e"+ 1 ) = 2At(Vp(t n+1 ),^ - e n ) < - e"|| 2 + 2(At) 2 \\Vp(t n+1 )\\ 2 . (37) 

Combining the inequalities (fM]) . (pJ5"j) . (|56")) . (|57|) . and rearranging terms, we obtain 

Ik™^ 1 !! 2 - l|e"|| 2 + ^Ik"^ 1 - e"|| 2 + ^Ai||Ve^+ 1 || 2 

< 2(Ai) 2 ||Vp(Wi)|| 2 + CAt\\e n \\ 2 + C(At) 2 (J*" +1 t\\u tt \Wdt + \\u t \\ 2 dt). (38) 

The step 4 of the algorithm (|T9" ]) -([2"T ]) yields 

e" = (1 - X )e" + X F(w n+1 )e n + x («(t n ) - F(w n+1 )u(t n )). (39) 

The definition of the filter and recalling that e™ is the L 2 projection of e n give j|F(w™ +1 )e™|| < ||e n || < ||e™||. 
We use this to deduce from ([59")) the following estimate: 

||el = (l-x)||el+ X ||F(^Vll+xN*r.) - F(w n+1 )u(t n )\\ < \\e"\\+x\Ht n )-F(w n+1 )u(t n )\\. 

Now we apply (|28[1 and square the resulting inequality to get (for the sake of convenience we assume At < C 
and recall x = XoAt): 

||e"|| 2 < (1 + A*)||e-|| 2 + CAtC ax . (40) 
We substitute HDl to the left-hand side of d38j for lie" II, use lie" II < ||e"|| and arrive at 



| e «+l||2 _ || e n||2 + || £ n+l _ £ , l+ l||2 + l|| e n+l _ e "|| 2 + uAt\\ Ve"+! || 2 

< 2(Ai) 2 ||Vp(t„ +1 )|| 2 + CAt\\^\\ 2 + C(At) 2 (^J^ tWuuf^dt + \\u t \\ 2 dt 
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Summing up (|4"Tj) from n — to n = I — 1, assuming that w° — w = uq (this implies ||e°|| = ||e°|| = 0), we 
obtain 

\w\\ 2 + J2 ii e " +1 - e " +1 " 2 + 5 Yl ii e " +1 - e "ii 2 + J2 ^ A iii Ve " +i n 2 

n— n— n— 

< CAt\\^\\ 2 + 2(Ai) 2 J2 II W„+i)l| 2 + C(Ai) 2 ( r t\\u tt \Wdt + T |Kfdi) + C5l 

n=0 n=0 ^ to 



S 2 

max 



i-1 



< ^ CAt||e"|| 2 + CAt + C5 2 m: 



n=0 



Applying the discrete Gronwall inequality yields (ptTj) . ID 

Now, we will use the result of the lemma and improve the predicted order of convergence for the velocity. 
The main result in this section is the following theorem, stating that all w n+1 , w n+1 and u n+1 are first-order 
approximations to the Navier-Stokes solution. 

Theorem 3 Assume the solution to the Navier-Stokes system satisfies |25|) and \ = XoAt. Suppose dQ. € C ' 
or is convex. It holds 

At ^(|k"|| 2 + ni 2 + ini 2 ) < C((At) 2 + CJ- (42) 

n=l 

Additionally assume ||Vpt|| 2 < C and the filtering radius is bounded as ^ ax < C At, then p n is an approx- 
imation to p(t n ) in L 2 (0)/R in the following sense: 

i 

AtJ2\\p n -p(tn)\\ 2 <C(At + S 2 niax ). (43) 

n=l 

Proof. Literally reaping the arguments from [22], pp. 66-69, one shows the estimate 

lk" +1 H 2 v - l|e n ||^ + - e n \\ 2 v , + ^At|| e n+1 || 2 < c(Ai||e" +1 || 2 y , 

+ (At) 2 J*" +1 (tWualW + \\u t \\ 2 )dt + (At) 2 \\Ve^\\ 2 + At\\e^ - e n \\ 2 + At\\e n+1 - e^f) . (44) 

The estimate ([29]) gives ||.Fe n ||y' < ||e"||y + <5 2 iaxll^ 7e "II- Here and in the rest of the proof the filtering is 
based on the w n+1 velocity, that is F- := F(w n+1 )-. Due to the assumption 30 <E C 1 ' 1 or O is convex, the L 2 
projection on H is H 1 stable, i.e. ||Ve n || < C||Ve™|| and therefore we conclude 

\\Fe n \\ v , <\\e n \\ v , +CS 2 max \\V^l 
Using this and (|2"9")1 , we get from (|31?|) for x = XoAt 

l|e n ||y< = (1 - xWWv + x\\Fe n \\v> + x\\u(t n ) - Fu(t n )\\ v , < \\e n \\ v , + CAt (<£aJ|Ve«|| + \\u(t n ) - Fu(t n )\\v) 

< ||e"|| v ,+CAt^ ax (||Ven|| + l). 
Squaring the inequality, we get after elementary calculations 

\\e n \\ 2 v , < (1 + At)\\e n \\ 2 v , + CAtSi ax (||V £ »|| 2 + l) . 
We substitute the above estimate to the left-hand side of (|44|) and arrive at 

\\e n+1 \\v< - lk"llv' + V l+1 - e n fv + ^Ai||6" +1 || 2 



< c(At(\\e n+1 \\ v , + \\e n f v ,) + (At) 2 jf " (t\\u tt \\\ + \\u t f)dt + (At) 2 ||V e »+if 

+ Atdl^+i - e "|| 2 + ||e" +1 - e^i|| 2 ) + Afe&Jl + l|V £ "|| 2 ) S 
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Assume for the sake of convenience 5 max < C. Summing up the inequalities for n = 0, . . . , I — 1, we get 

l-i l-i 



ll^ + Ell^-^llv' + E"^!!^ 1 !! 2 

< C {y, ^¥ n+1 \\v> + (At) 2 f\\\uu\\ 2 v . + IM^ + ^E A*l|Ve«|| 2 

\n=0 n=0 

l-l l-l \ 

+ ^Ai|| e ™+i- e ' l || 2 + E At ll £ ' l+1 - e ' l+1 ll 2 + A^ ax • (45) 

n=0 n=0 / 

Now we use the result of the Lemma [3] to bound 

l-i i-i i-i 

+ C E At||V6«+i|| 2 + E At||e" +1 - e"|| 2 + £ At|| e " +1 - e«+i|| 2 < C((At) 2 + Ai^ + O- 

n=0 n=0 n=0 

Thus, applying the Gronwall inequality to (|45p yields 

i-i i-i 
He'll?,, + E H £ " +1 " e "Hv' + E ^|| e " +1 || 2 < C((At) 2 + O- (46) 

n=0 n=0 

Here we also used AiJ 2 ^ < (At) 2 + ^ ax . Finally, the Lemma [3] helps us to estimate 

l-i i-i i-i 



A* ^ n e «+iir < At j2 - t n+i \\ 2 + At e ii e " +1 n 2 ^ c ((a*) 2 + O- 

n— n— n— 

/ i-1 Z-l 

At^ ||e"|| 2 < At£ - e"|| 2 + At£ ||e" +1 || 2 < C((At) 2 + < ax ). 



n—Q n—Q n—Q 

These estimates together with proves the velocity error estimate of the theorem. 

Further we show that the pressure is weakly | order convergent to the true solution. Denote the pressure 
error as q n = p n — p(t n ). We may assume (q n , 1) = 0. It holds 

- Vq n+1 = -^ n+1 - e n ) + ^Ae^+i + (w n ■ V)^ 1 - (u(t n+1 ) ■ V)u(t n+1 ) + R n . (47) 

Repeating the arguments from [32] and using the Necas inequality, see [37J , one deduces from (jT7|) 

(Vn n+1 v\ 1 ■ 

|| g » +1 || < c sup [ q - ' > < — 1| £ » +1 - e»||_ a + C(\\R n U + ||Ve«+i|| + ||V e » +1 || + \\u(t n+1 ) - u(t n )\\). 
veH^nr l|Vu|| At 

Therefore, by using (j3Tj) . we get 

i-i i-i 

At \W l+1 II 2 < ^ E H V ( £ ' l+1 - Oll-l + C ( At + ^ (48) 
n=0 n=0 

To bound the first term on the right-hand side of (|48[) one estimates: 

|| e «+l _ e n||_ i < c || e n+l _ e n|| < c (|| e «+l _ e ™|| + || e « _ e "||) < c (||^+i - e »|| + || e « _ e "||). (49) 

The estimate for the second term on the right-hand side of (|49|) follows from (|39|) : 

||e™ - e"|| < XoAt(||e" - Fe"|| + \\u(t n ) - Fu(t n )\\) < XoAt(||e n || + \\Fe n \\ + \\u(t n ) - Fu(t n )\\). 
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Thanks to (|28j). (|3T|) . and H-Fe"!! < ||e ra || we continue the above estimate as 

||e"-e n || <C((At)i + At<5 max ). (50) 

Below we shall prove the bound 



l-i 



J2h n+1 -t n \\ 2 <C((Aty + At6i ax ). 



n=0 



From dHJ) and ([21]) we get 



i_( e n+l _ e «) _ + V p(t„ +1 ) + (W n ■ V)^ 1 - (U(t n+1 ) ■ V)u(t n+1 ) = W\ (51) 

The projection step (gDJ) gives e n = e" + AtVp n , so dHH]) yields 

e" = (1 - + AiVp") + X Fe n + X (u(i») - Ftt(t„)). 

Substituting this in (|5ip implies 

i_ (e ^i _ e n) _ ^Ae^i + (1 - X )V(p(i„ +1 ) - p ») + x Vp(i n+1 ) - |^(Fe" - e") - ^(u(t n ) - Fu(t n )) 

+ (w n ■ V)w^ 1 - (u(t n+1 ) ■ VMVn) = J? n . (52) 

The inner product of ([52"]) with At(e n+1 — e n ) gives 

_ e n||2 + ^(|| Ve ^l||2 _ || Ve n||2 + || V (^+i - e ™)|| 2 ) 

= At(ir, e^ 1 -e«) + (1 -x) At(p(t n+1 ) -p n , div(e^ -e»)) + At((u>" • V)^ 1 - ( M (t Il+1 ) • V) U (t„ +1 ), e^ 1 -e^) 

- xA<(Vp(t„+i), e^+i - e") + x(^e™ - e™, e™^ 1 ™ e Tl ) + x(u(«n) - Fu(t„), e"^ 1 - e") 
= At(J? n ,e^ - e») + (1 - X )At [(<?", div^ 1 - e")) + (p(*„+i) - p(t„), div^ 1 - e^)) 

- x [At(Vp(t n+1 ), e^ 1 - e") - (Fe n - e", e^ 1 - e") - (u(t„) - Fit(t n ), e^ 1 - e") 

+ At((w n • V)^ 1 - («(tn+i) ■ V)«(tn+l), - e") 

= 7 1 +/ 2 +l3+74+/5+-f6+^7- (53) 

The last term J7 is estimated in [26] : 

At|((io" • V)^ 1 - («(t„ + i) • V)u(t n+1 ),^ - e*)| 
< cr||^ i+1 — ^|| 2 + C((ZXt) 2 1| 2 + (ZXt) 2 ||e" +1 || 2 + ||Ve"|| 2 ||V?" +1 || 2 + ^^ 1 ||V(e^+ 1 — ^)l| 2 + (^*) 3 ) 

for some a > 0, which can be taken sufficiently small. Applying ([3"T]) and ||Ve"|| < C||Ve n || leads to 

h < oW +1 ? l f + C((Atf + (At) 2 8 2 max ) + (At)l ||Ve»]| 2 ||W +1 || 2 + ^||V(e^ - e»)|| 2 . (54) 
For J4, 1$, and Iq one has 

h = - X At(V P (t n+1 ),e^-^) < C X 2 {At) 2 \\V P (t n+1 )\\ 2 + a\\^+ 1 -e~"|| 2 , (55) 
4 = X(^e" - e",^ 1 - f~ n ) < C X 2 (\\Fe n \\ 2 + \\^\\ 2 ) + a\\^ - e^|| 2 

< C((At) 3 + (A^LJ + CT ||e^ - £ »|| 2 , (56) 

h = X(u(t n ) - Fu(t n ), - e") < C(At) 2 ^ ax + - e~>|| 2 . (57) 
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The terms h, I2 and I3 are estimated in [22]. Using those estimates and (|54| -(|57 j) in (|53|) yields for sufficiently 
small a > 0: 

||^+i _ e n||2 + ^(nve^f - ||Ve^|| 2 ) + (1 - ^(A^dlV^ 1 !! 2 - j|V<f || 2 ) 

< C {(At) 2 f n+1 \\u tt \\ 2 dt + (At) 2 [ tn+1 \\V Pt \\ 2 dt + (At) 4 \\V P (t n+1 )\\ 2 



+(Atf + (A^Lx + Atl||Ve«|| 2 || W l+1 || 2 } . (58) 



We sum up the estimate for n = 0, . . . , I — 1 and apply our assumptions for the solution to Navier-Stokes 
solution. This leads to the bound 

£ ll^+i _ e ~«|| 2 + ^|| V ?|| 2 < C((At) 2 + Ai<£ ax + (At)* ^ ||V e ~i 2 ||V^ +1 || 2 ). 

n=Q n=0 

The application of the discrete Gronwall inequality, (|31l) and the assumption <5^ lax < CAt yields 
£ ll^+i _ e n|| 2 + ^||v?|| 2 < C ((Ai) 2 + A^Lx) exp {(At)* f] ||Ve^f] 

n=0 I n=0 J 

< C ((At) 2 + At^ax) cx P {ff((At)* + (At)-»0} 

< C((At) 2 + A*0- 
Therefore, (|4"8 ]l -([50 |) yield the desired bound: 

At-£\\q n +Y<C(At + S 2 max ). 

n=Q 

u 
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